Memory effects in non-adiabatic molecular dynamics at metal surfaces 
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We study the effect of temporal correlation in a Langevin equation describing non-adiabatic 
dynamics at metal surfaces. For a harmonic oscillator the Langevin equation preserves the quantum 
dynamics exactly and it is demonstrated that memory effects are needed in order to conserve the 
ground state energy of the oscillator. We then compare the result of Langevin dynamics in a 
harmonic potential with a perturbative master equation approach and show that the Langevin 
equation gives a better description in the non-perturbative range of high temperatures and large 
friction. Unlike the master equation, this approach is readily extended to anharmonic potentials. 
Using density functional theory we calculate representative Langevin trajectories for associative 
desorption of N2 from Ru(OOOl) and find that memory effects lowers the dissipation of energy. 
Finally, we propose an ab-initio scheme to calculate the temporal correlation function and dynamical 
friction within density functional theory. 
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I. INTRODUCTION 

Modern computational surface chemistry, as e.g. ap- 
plied to heterogeneous catalysis, is largely based on the 
Born-Oppenheimer approximation and potential energy 
surfaces which are typically obtained using density func- 
tional theory (DFT).[2] In the adiabatic approximation 
the electrons are assumed to follow the motion of the nu- 
clei instantaneously and the dynamics thus becomes con- 
fined to the ground state potential energy surface. While 
the adiabatic approximation has certainly been success- 
ful in giving a detailed quantitative account of a range of 
chemical reactions on metal surfaces, it is still not clear 
under which general circumstances the approximation is 
reliable. In particular, the role of non-adiabatic ef- 
fects is often difficult to asses due to the inadequacy of 
low dimensional models of surface dynamics. For exam- 
ple, unusual sticking coefficients in the measured disso- 
ciative adsorption of N2 on Ru(OOOl) hints at strong 
non-adiabatic energy loss, 6] but has been accounted 
for by multi-dimensional adiabatic dynamics. 0, Q For 
other reactions, such as associative desorption of N2 
from Ru(OOOl), non-adiabatic effects still seem to be very 
important |[2|, |9|, ilfl] and multi-dimensional adiabatic sim- 
ulations have not been able to account for large energy 
losses during desorption. [il| Another example where adi- 
abatic dynamics have failed is the dissociation of O2 on 
Al(lll) where spin selection rules gives rise to highly 
non-adiabatic behavior. (T2j 

Non-adiabatic dynamics for isolated molecules is usu- 
ally handled by including the first few excited adiabatic 
potential energy surfaces and imposing some surface hop- 
ping scheme. When distinct adsorbate diabatic states 
are present there may be physical arguments why the 
adsorbate should remain in such a state during a reac- 
tion and the non-adiabatic dynamics can be evaluated 
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by constraining the adsorbate to such a diabat.[13] This 
picture may then be improved by introducing surface 
hopping between diabats. However, for molecules ad- 
sorbed on metal surfaces there is an infinity of electronic 
excited states in the immediate vicinity of the ground 
state and surface hopping may not be the most practi- 
cal scheme. Another popular and rather general method 
to handle non-adiabatic effects is through Langevin dy- 
namics where electronic friction and stochastic forces ac- 
count for dissipation and fluctuation as a result of cou- 
pling to excited electronic states. [l^ - [l7j Usually, the so- 
called Markov approximation is employed where the fluc- 
tuating forces are not temporally correlated but can be 
related to the electronic friction through the fluctuation- 
dissipation theorem. [ISj] At high electronic temperatures, 
thermal excitations dominate the electronic system and 
the Markov approximation is good for describing chemi- 
cal reactions mediated by hot electrons, [l 9, 20] However, 
for non-adiabatic dynamics in general, the Markov ap- 
proximation may fail and it then becomes important to 
take into account the 'memory' of the system. 

In the present paper we explore the consequences of 
Langevin dynamics with and without the Markov ap- 
proximation at various tem per atures. We follow the ap- 
proach of Brandbyge et al. [I4I and base the analysis on 
a model Hamiltonian from which the electronic friction 
and correlation function can be derived explicitly. We 
start by modelling the internal stretch mode of CO ad- 
sorbed on Cu(lOO) by a considering a harmonic oscillator 
coupled to a thermal reservoir of electrons and compare 
the results to those obtained with a master equation ap- 
proach. The general trend we see is that at low temper- 
atures the Markov approximation overestimates the ef- 
fect of dissipation. This is because the fluctuating forces 
are of thermal origin in the Markov approximation while 
the dissipative terms originate from non-thermal excita- 
tions and the relative effect of dissipation compared to 
fluctuations is thus increased. We also show that non- 
Markovian dynamics is needed in order to ensure energy 
conservation at low temperatures and thus maintains de- 
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tailed balance between fluctuations and dissipation. As- 
sociative desorption of N2 from Ru(OOOl) is then studied 
using the Langevin equation on representative trajecto- 
ries and again, memory effects are shown to reduce the 
dissipation of energy. Finally, we comment on a possible 
method to obtain the full correlation function and thus 
include memory effects in an ab-initio setting based on 
DFT.^21.] In appendix \X[ we review the connection be- 
tween the reduced density matrix and the Langevin equa- 
tion and emphasize the probabilistic interpretation of the 
correlation function. In appendix |B1 we review how cor- 
related stochastic forces can be randomly sampled given 
a discretized version of the correlation function. 



the resonant state \a). However, usually we are only in- 
terested in the nuclear degrees of freedom and it is then 
convenient to trace out the electronic degrees of freedom 
from the full dynamics. This is accomplished by the re- 
duced time dependent density matrix: 

p,e<iW = Trei(e-*^*/Voe^^*/''), (4) 

where Tt^i means the trace over electronic states and po is 
the full density matrix at t — to- Choosing an adsorbate 
basis the diagonal elements of the reduced density 
matrix give the probabilities that the adsorbate is in a 
particular state at time t. 



II. MODEL 



A. Non-Markovian master equation 



A commonly used electronic Hamiltonian describing 
of atoms or molecules adsorbed on metals surfaces is the 
Newns- Anderson model, 0, [2^ where the adsorbate is 
described by a single adsorbate state \a) which hybridizes 
with metallic states \k) and thus acquires a broadening 
in energy. A very simple non-adiabatic extension of this 
model is obtained by coupling the resonant states \a) 
to an adsorbate degree of freedom x and extending the 
Hamiltonian with a nuclear kinetic energy and adiabatic 
potential. Assuming a quadratic nuclear potential and 
linear coupling to the resonance, the Hamiltonian be- 
comes 



H — Hei + i?o + Hi, 



(1) 



k k 

Hi = -fclcax, 

where p is the nuclear momentum, M the adsorbate ef- 
fective mass, and and cj. are creation operators for ad- 
sorbate and metallic electronic states respectively. The 
Hamiltonian Hq + Hi can be thought of as a harmonic 
oscillator which is shifted when the resonance becomes 
occupied and the coupling constant / is the force felt 
by adsorbate in this state. We will furthermore restrict 
ourselves to the wide band approximation in which the 
metallic band of electrons is assumed to be much wider 
than the width of the resonant state. The resonance pro- 
jected density of states is then a Lorentzian: 
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7r(e-eo)2 + (r/2)2' 
with full width at half maximum given by 



r = 27r 



k 



(2) 



(3) 



Due to the non-adiabatic coupling in Eq. [U the adsor- 
bate may exchange energy with the electronic system via 



We will consider the time-dependent probability of be- 
ing in a particular energy eigenstate |n). The equation 
governing these probabilities is known as a master equa- 
tion and can be derived by taking the trace of the Li- 
ouville equation for the full density matrix. The result 



dpr 



dt 



-[Ho,Pr 



(5) 



where the influence functional = TTei[Hi, p]/h de- 
pends on the complete history of the full density ma- 
trix. Gao [13] has shown how to evaluate J-[p] using the 
Hamiltonian ([T]) within the self-consistent Born approxi- 
mation. Furthermore, imposing the Markov approxima- 
tion, where it is assumed that the influence functional 
only depends on the instantaneous value of the density 
matrix, and taking the diagonal elements of Eq. ([S]) led to 
an explicit expression for the master equation. However, 
using the formalism of Gao ^23], it is straightforward to 
generalize the results to a non-Markovian master equa- 
tion. For completeness we state the result here which 
is 



dt' 



Wn^^{t-t')pn{t') (6) 



Wm^n{t-t')p^{t') , 



with the differential rates given by 

Wn^m{t) = / doJi / (iW2Pa(w2)(l - «f(w2)) (7) 



X pa{uJi)nFiuJl) C0S[(W1 -UJ2 + LJnm)t], 

where Pmit) — {m\pred(t)\m), \m) is an eigenstate of Hq 
with eigenvalue Xmn — (m|a;|n), np{e) is the Fcrmi- 
Dirac distribution, pa{£) is the projected density of states 

^ and UJnm = (Sn - £m)/h- 

The Markov approximation is obtained by extending 
the temporal integration to infinity which is justified 
when t is much larger than some electronic correlation 
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FIG. 1: The differential transition rate Wo-)-i(i) given in 
Eq. (O using three different temperatures and F = 2.0 eV , 
£0 = 2.6 eV, hwo = 0.250 eV and / = 8.7 eV/k. The figure 
shows that the time dependence vanishes after a few fem- 
toseconds and since the typical timescale of change in Pnit) 
is ~ 100 /s, the Markov approximation is expected to work 
well for the master equation. 



time tc- The probabilities Pm are then assumed to de- 
pend on t rather than t' and integrating over t' yields 
the usual golden rule type expression for the transition 
rates. [24j The master equation was derived assuming that 
ofF-diagonal elements of the density matrix are not im- 
portant. Allthough it is straightforward to generalize the 
result (|n]) to a coherent master equation which takes off- 
diagonal elements into account, it has been shown that, 
if the initial state is not coherent, the off-diagonal el- 
ements will have very little influence on the diagonal 
elements. [2^ In terms of multiple inelastic scattering 
events, neglecting coherency corresponds to associating a 
probability distribution to each scattering event. [IH US] 
In Fig. [1] we show the differential transition rate 
Wo-i.i(i) at three different temperatures. The time de- 
pendence only depends on the properties of the electronic 
system at the given temperature and the non-adiabatic 
coupling / simply gives an overall scaling. Since the prob- 
abilities Pn{t) typically change on timescales of ^ 100 fs 
and the Wn^m{t) approach zero within a few femtosec- 
onds, the Markov approximation is expected to be very 
good for the master equation in a large range of temper- 
atures. 



B. Non-Markovian Langevin dynamics 

If we calculate the diagonal of the reduced density ma- 
trix in a basis of position eigenstates, a Langevin equa- 
tion emerges. This is achieved by writing Eq. as a 
path integral and using the Feynman- Vernon formalism 
of influence functionals to obtain an effective action to 
second order in the frictional coupling [13 The 

result is given in Eq. (jA6l) . This approach is not pertur- 
bative in the same sense as the master equation where 



the derivation is based on a direct second order expan- 
sion of the reduced density matrix. Rather, the second 
order expansion of the action leads to a density matrix 
which contains all orders of the frictional coupling. As 
explained in appendix [Xj the result can be interpreted 
as a sum over classical Langevin trajectories with initial 
conditions sampled from the Wigner distribution of the 
initial state and the equation of motion is 



Miit) 



Muj'^x{t) 



dt'r]{t-t')x{t')=i{t), (8) 



where •q{t) the dynamical electronic friction and ^(t) is 
a Gaussian distributed stochastic force specified by its 
correlation function: 



mi{t')) = K{t-t'). 



(9) 



With the model Hamiltonian ([T|) it is possible to eval- 
uate the friction and correlation function explicitly. The 
result is: [13 



2^ 



A(w) cos{ijjt), 



(10) 



with 



A{0J) = - dw2G{u,,U2) (11) 

^ J -00 ^'^ J-00 

X (5^0; - (w2 - uji)^ (npiuJi) - ?iF(w2)), 



G{UJI,UJ2) = 4:71 f Pa{i^l)Pa{i^2) 



(12) 



and Pa{i^) is the projected density of states Eq. ([2]). The 
correlation function is 



Kit) 



n 



duj 



■a;A(w)cothf ^ cos(a;i). (13) 

27r ^ ' \2kBTJ ^ ' ^ ' 



In Fig. [2] we show the correlation function for three 
different temperature. The structure and typical correla- 
tion time is very similar to the differential rate shown in 
Fig. [T] However, in contrast to the master equation the 
timescale of motion in the Langevin equation is on the 
order of At ^ 1 fs and correlation effects may be very 
important at low temperatures. In general the stochastic 
forces at a given time depends on the state of the elec- 
tronic system which again depends on the path taken 
by the adsorbate. This gives rise to correlation between 
forces at different times and this "memory" is the price 
one pays for tracing out the electronic degrees of freedom. 
The range of memory in the system depends strongly on 
temperature since high temperatures tend to rapidly de- 
stroy coherence in the state of the electronic system. In 
the high temperature limit where fc^T ^ h/ At, one ob- 
tains the well known expression 



K{t) = 2kBTf]oS{t), 



(14) 
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FIG. 2: The correlation function K{t) given in Eq. (|13|l using 
three different temperatures with F = 2.0 eV, eo ~ 2.6 eV, 
hujo = 0.250 eV and / = 8.7 eV/K. At low temperatures 
correlation can persists for several femtoseconds. 



FIG. 3: The correlation time Eq. (jlSp as a function of tem- 
perature with. Below 3500 K the correlation time becomes 
larger than2 a femtosecond which is the largest time step we 
can use in the molecular dynamics and non-Markovian pro- 
cesses therefore begins to influence the dynamics below this 
temperature. 



where rjo = A(0)/2. This is the Markov approximation in 
which there is no correlation between forces at different 
times. Taking CO on Cu(lOO) as an example, the small- 
est timescale is the period of vibrational motion which is 
^ 16 fs. With a standard Verlet integration one needs a 
timestep of At ~ 1/s to describe ground state vibrations 
and a first estimate of the validity of the Markov approx- 
imation is obtained as: T > h/{AtkB) 2900 K. To 
get a better quantitative estimate of the validity of the 
Markov approximation we can consider the correlation 
time tc given by 



JdtK{t) 



(15) 



It should be noted that the correlation time is only a 
function of the electronic system and does not depend 
on the non-adiabatic coupling /. We have calculated 
tc as a function of temperature and the result is shown 
in Fig. [31 For molecular dynamics requiring a time 
step no larger than ~ 1 fs, we see that the correla- 
tion time becomes larger than this when the tempera- 
ture comes below 3500 K. Thus below this temperature 
non-Markovian processes play an important role in the 
dynamics. 



III. RESULTS 

Before we test the role of non-Markovian effects on a 
generic non-adiabatic surface reaction, we will compare 
Markovian and non-Markovian Langevin dynamics for a 
harmonic oscillator potential with results obtained from 
a master equations approach. 



A. Quadratic potential 

It is easy to see that the fluctuating force in the 
Langevin equation has to vanish within the Markov ap- 
proximation Eq. (IT4l) when T^i — > 0. With a quadratic 
potential and r]{t) = A{0)/26{t) it is then possible to 
solve the Langevin equation analytically which gives the 
time-dependent energy 



EMarkovit) = -Boe"*/", T = 2M/A(0), 



(16) 



where Eq is the initial energy. However, as shown in Ref. 
[20l . the Langevin equation is quantum mechanically ex- 
act for a harmonic potential if the initial conditions are 
accounted for correctly and the total energy should thus 
not be allowed to decay below huj/2. The problem is 
that the Markov approximation neglects all non-thermal 
excitations of the electronic system and leads to pure 
dissipation at T^i — > 0. In reality, an oscillating adsor- 
bate will induce (non-thermal) excitations of the elec- 
tron gas which may then influence the propagation of 
the adsorbate. In general, it is therefore expected that 
the Markov approximation tends to underestimate the 
influence of the electronic system on the adsorbate. This 
non-Markovian effect should vanish at high temperature 
where the thermal excitations of the electronic system 
dominate. In Fig. U we show the time evolution of the 
average energy of a harmonic oscillator interaction with 
a thermal reservoir of electrons at six different temper- 
atures. The average energy is calculated using the full 
non-Markovian correlation function as described in ap- 
pendix |B] and within the Markov approximation. The 
initial state was chosen as the vibrational ground state 
and included exactly by phase space sampling the Wigner 
distribution. [2^ The parameters used were chosen to 
match the internal vibrational mode of CO adsorbed on 
Cu(100)[23, ill and we have thus taken F = 2.0 eV, 
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FIG. 4: Average energy of a harmonic oscillator interacting 
with a thermal reservoir of electrons at six different temper- 
atures evaluated using Langevin dynamics with memory and 
with the Markov approximation. The Markov approximation 
fails below T = 3000 K where quantum fluctuations are im- 
portant. 



eo = 2.6 eV, huj ^ 0.25 eV, and / = -8.7 eV/A. The 
failure of the Markov approximation and resulting decay 
of the average energy is clearly seen at low temperatures. 
In particular, at T = 500 K the Markov approximation 
gives rise to exponentially decaying energy whereas the 
energy remains nearly fixed dX E k, when memory ef- 
fects are included. For high temperatures (T > 3000 K) 
thermal excitations dominate and the Markov approxi- 
mation becomes reliable. In all calculations we have con- 
verged the results by decreasing the time steps. 

In Fig. [5] we show the average energy of the harmonic 
oscillator after 5 ps of interaction with a thermal reser- 
voir of electrons. The energy is calculated with non- 
Mar kovian Langevin dynamics, Mar kovian Langevin dy- 
namics, the master equation with rates obtained from 
perturbation theory and the master equation with ex- 
act (non-perturbative) rates, [l^] Using the full non- 
Markovian master equation Eqs. ©-(O does not change 
the results. The non-Markovian Langevin approach 
matches the exact non-perturbative Master equation ap- 
proach, whereas the perturbative Master equation fails 
at high temperatures and the Markovian Langevin ap- 
proach fails at low temperatures. It may be surprising 
that the non-Markovian Langevin equation reproduces 
the exact and not the perturbative master equation. 
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FIG. 5: Average energy of a harmonic oscillator after 5 ps 
of interaction with a thermal reservoir of hot electrons. The 
non-Markovian Langevin equation and the non-perturbative 
master equation both give the correct dependence, whereas 
the Markovian Langevin equation fails at low temperatures 
and the perturbative master equation fails at high tempera- 
tures. 



However, as mentioned above, the perturbative deriva- 
tion of the master equation is based on a direct evalua- 
tion of the reduced density matrix to second order in the 
non-adiabatic coupling, 24'] while the Langevin equation 
is derived by constructing an effective action to second 
order in the non-adiabatic coupling. Thus, while the 
reduced density matrix calculated from the effective ac- 
tion only becomes exact in the small friction limit, it does 
contain terms to all orders in the non-adiabatic coupling 
and is a much better approximation for large frictional 
coupling and high temperatures than the direct pertur- 
bative derivation leading to the master equation Eqs. 
0. 

It may seem like a complete overkill to apply a non- 
adiabatic Langevin dynamics to a harmonic potential 
when the results are readily obtainable from the master 
equation approach. However, for anharmonic potentials 
it is not possible to derive transition rates for the master 
equation exactly and the best approximation is then the 
non-Markovian Langevin equation. This was also con- 
cluded in Ref. [l^ where a perturbative master equation 
approach underestimated transition rates in a Morse po- 
tential compared to a Markovian Langevin approach. 



B. Associative desorption of N2 from Ru(OOOl) 

As an example of a potential where the master equa- 
tion approach is not readily applicable, we consider the 
well-known example of associative desorption of N2 from 
Ru(OOOl). This system has been subject to extensive 
experimentally, [^, and theoreticaljl 0-11 HH studies 
and much evidence points to a non-adiabatic dissipation 
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of energy during associative desorption. 

The Langevin equation can be generalized to an ar- 
bitrary potential Vo{x) by a semiclassical expansion of 
the potential and the excited state forces acting on 
the adsorbate. J7j The potential is included by making 
the substitution Muiqx'^/2 — s> Vo{x) in the Hamiltonian 
Eq. ([1} and the friction arises from an excited resonant 
state with potential energy Vi{x) and is included by 
the generalizations —fx -> Sa{x) = Vi{x) — Vo{x) and 
Vak ~^ Vakix) which in the wide band limit leads to a 
position dependent resonance width r(a;). When mul- 
tiple coordinates {xt} are considered 77, A, and K{t) in 
Eqs. (|T0l) - (fT3)) become tensors and the Langevin equa- 
tions for each coordinate become coupled through terms 
like rjij {x)xj . In addition to temporal correlation, the 
stochastic forces acting on different coordinates become 
correlated through the off-diagonal terms in the corre- 
lated function: 

umm^K^At)- (17) 

In fact, since the friction tensor is well approximated by 
Aij{uj) (X fifj with fi = d£a{x)/dxi, Kij(t) has only 
one non-zero eigenvalue. This implies that there is a sin- 
gle (coordinate dependent) mode on which the stochastic 
force acts and the random forces can thus be regarded as 
completely spatially correlated at any given time. 

We have studied associative desorption of N2 from 
Ru(OOOl) using the code gpaw,[2^ [30] which is a real- 
space Density Functional Theory (DFT) code that uses 
the projector augmented wave method. [Sll, [s^ The 
Ru(OOOl) substrate was modelled by a three layer slab 
where the top layer was relaxed. We used a grid spacing 
of 0.2 A and the calculations were performed in a (2x4) 
supercell sampled by a (4x6) grid of fc-points using the 
RPBEfs^ exchange-correlation functional. The friction 
is assumed to be dominated by the 27r orbital which is 
only partly occupied in the ground state. To calculate the 
excited state potential energy Vi{x) we applied a gener- 
alization of the A-self-consistent field method where the 
resonant state is expanded in a basis of Kohn-Sham or- 
bitals and the resulting resonant density is added to the 
density in each iteration step. Thus for each adsorbate 
position we calculate the energy resulting from forcing 
an electron into a 2tt orbital. For details on the method 
and comparison with experiments we refer to Ref. [s^ 
We have restricted the analysis to the two-dimensional 
desorption process considered in Ref. where the two N 
atoms are adsorbed at adjacent hep hollow sites and des- 
orbs by moving perpendicular to the bridge towards the 
fee hollow while changing the center of mass coordinate. 
While a two-dimensional analysis is almost certainly not 
sufficient to obtain quantitative results, 0, [Tlj we do 
expect to draw some qualitative conclusions about the 
validity of the Markov approximation for this system. 
The calculated ground and excited state potential en- 
ergy surfaces are shown in Fig. [5] To obtain r((i, z) we 
have fitted the width of the projected density of states of 
the 27r orbital along the minimum reaction to an expo- 
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FIG. 6: Ground and excited state potential energy surfaces 
of N2 adsorbed on Ru(OOOl). The excited state was obtained 
by occupying the 2-k orbital of N2. 

nential r(z) = Foe-^^"^')/^'* and obtained Fq = 3.0 eV, 
Zd = 0.5 A and zt is the center of mass position at the 
transition state. The frictional force dea{d,z)/dd in the 
internal mode become large in the exit channel and gives 
rise to large dissipation of the internal energy while the 
molecule desorbs. However due to the rapid decay of F(z) 
the friction tensor essentially vanishes at 2 = 3 A. The 
amount of dissipated energy thus largely depends on the 
time spend in the exit channel in the immediate vicinity 
of the transition state. 

To examine the impact of non-adiabatic dissipation of 
energy and, in particular, the validity of the Markov ap- 
proximation, we have considered four representative ini- 
tial conditions leading to desorption. All four are initially 
at the transition state with a kinetic energy of 0.1 eV. 
The kinetic energy is then concentrated in positive or 
negative center of mass momentum or positive or nega- 
tive internal momentum. We have taken the surface and 
thus the electronic temperature to be T = 900 ii'.fl], [§]. 
Table |I] displays the average energy loss in a desorp- 
tion event of the four initial conditions with and with- 
out the Markov approximation. The reason for the large 
difference is due to the average time spend in the exit 
channel which for initial negative internal momentum is 
^ 125 fs and for initial positive center of mass momen- 
tum is ^ 250 fs. The shift to lower dissipation in non- 
Markovian dynamics is what we would expect from the 
conclusions in Sec. IIII Al and Fig S) In general, memory 
effects tends to increase the importance of fluctuating 
forces and thus decrease the overall dissipation of energy. 

The present analysis should in no way be taken as 
a quantitative study of non-adiabatic effect in associa- 
tive desorption of N2 from Ru(OOOl). In such a study 
one would need to sample a thermal distribution of ini- 
tial configurations at the transition state and include 
all 6 molecular degrees of freedom. [ll| Furthermore, the 
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Mode 


z- 


z+ 


d- 


d+ 


Markovian 


0.31 


0.49 


0.079 


0.10 


non-Markovian 


0.13 


0.42 


0.053 


0.10 



TABLE I: Average energy loss (all numbers in eV) of trajec- 
tories leading to desorption for four different initial conditions 
with and without the Markov approximation at T = 900 K. 
The initial conditions were all at the transition state with 
a kinetic energy of 0.1 eV . d and z denotes initial momen- 
tum in the internal vibrational mode and the center of mass 
mode respectively and - and -f denotes the sign of the initial 
momentum. In general, the Markov approximation tends to 
underestimate the effect of fluctuating forces which results in 
too much dissipation. 



present study assumes that the electronic friction origi- 
nates from a single resonance (27r) and it is well described 
by the wide band approximation. That this is not the 
case has already been established and a complete ab- 
initio scheme for the electronic friction is needed. Such 
a scheme based on DFT has already been suggested and 
put to use within the Markov approximation0, [l^ [HI , 
but need to be generalized slightly to take memory ef- 
fects into account. In section ITVl we will propose such a 
generalization. 

For anharmonic potentials where the friction tensor ac- 
quires a position dependence, non-Markovian Langevin 
dynamics is rather time consuming, since one has to cal- 
culate and diagonalize the correlation function in each 
time step. In the present simulation, the time required 
for a single time step in the dynamics was increased by a 
factor of 10"^ compared to Markovian dynamics, but the 
computational time is, however, vanishing compared to 
that required for a full DFT calculation at a given posi- 
tion. For N2 on Ru(OOOl) the memory effects are clearly 
seen but probably not important compared to neglecting 
four degrees of freedom. The calculated energy dissipa- 
tion is not large enough to account for the vibrational 
damping observed in Ref. [lol but it is very likely that 
inclusion of more degrees of freedom would result in a 
larger amount of time spend in the exit channel and thus 
a much larger dissipation of internal energy. 



IV. NON-MARKOVIAN FRICTION AND 
FLUCTUATING FORCES FROM DENSITY 
FUNCTIONAL THEORY 

Within linear response theory, it is possible to derive 
an expression for the electronic friction for a general non- 
adiabatic Hamiltonian if one assumes classical adsorbate 
motion. [21] The result depends on the response function 
of the electronic system as well as the derivative of the 
electron-vibron coupling with respect to adsorbate co- 
ordinates. Replacing the true response function with a 
Kohn-Sham response function and the coupling Hamil- 
tonian by a Kohn-Sham potential, give the result for the 



electronic friction 

r; = -7rfi^|(^,|^|V^,}|' (18) 

f dnpie) 
X / de — 2 "{^i ~ ~ 

where ipi are Kohn-Sham orbitals with eigenenergies 
Ei. The result is valid within the Markov approxima- 
tion, since the memory in the Kohn-Sham potential 
has been neglected. Generalizing this result to include 
non-Markovian dynamics would require a non-adiabatic 
exchange-correlation potential, which is presently out of 
reach. However, since the result is equivalent to that ob- 
tained within the reduced density matrix formalism and 
the Hamiltonian ([T]) we can impose a very simple gen- 
eralization which reduces to the adiabatic result in 
the Markov approximation and to (fT0|) -([T3 | in the case of 
non-interacting electrons. Indeed, it is easy to verify that 
([T8| reduces to the Markovian limit of dTOJ-dlSl) if Vks 
is replaced with the Hamiltonian ([T]) and one is led to a 
non-Markovian Langevin equation based on DFT which 
is given by dHJ-dn]) and dH]), but with 

X S{uJi - £jh)5{uj2 - £j/K). 

This result for the dynamical friction was also derived in 
Ref. Ull, however, the dominating non-Markovian effect 
is the correlation function (|13p which follows from the 
reduced density matrix formalism in conjunction with 
(Hi). 



V. SUMMARY 

From a fundamental point of view it is important to 
realize that the Langevin equation gives an exact descrip- 
tion of a harmonic oscillator interacting with a reservoir 
of electrons if the initial quantum state is taken correctly 
into account. [20[ However, it is easy to see that the fluc- 
tuating forces vanish at low temperatures in the Markov 
approximation, which then results in a exponentially de- 
caying energy of the oscillator. This of course contradicts 
the quantum description of the oscillator and the prob- 
lem can be traced to the Markov approximation which 
does not take non-thermal electronic excitations into ac- 
count. In Fig. m we have shown explicitly how memory 
effects 'saves' the quantum behavior of the oscillator and 
conserves the energy of the vibrational ground state at 
low electronic temperatures. 

Another way of handling dissipative systems is using 
the master equation. This approach is based on a pertur- 
bative calculation of the reduced density matrix in basis 
of energy eigenstates. While this approach is fast and 
intuitively appealing, the method breaks down at high 
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temperature or large friction due to the perturbative na- 
ture of the method. In contrast, the Langevin equation 
is based on an effective action Eq. (jA2| - (jA3| giving a 
non-perturbative flavor. Furthermore, the master equa- 
tion requires quantization of the potential energy surface 
and becomes impractical for complicated potentials with 
many bound states. 

As an example of such a potential we have consid- 
ered the associative desorption of N2 from Ru(OOOl). We 
have not tried to perform a quantitative two-dimensional 
study of this system as done by Luntz et al. , but rather 
examined the effect of temporal correlation in two rep- 
resentative trajectories. As expected, the effect is an in- 
creased significance of the fluctuating forces leading to 
lower dissipation when memory is included. While this 
is may be of qualitative interest, the effect of including 
all degrees of freedom and performing an ab-initio cal- 
culation of the friction tensor, would almost certainly 
lead to corrections which are quantitatively much more 
important. jlT] 

Finally, we have provided an expression for the corre- 
lation function within an ab-initio DFT scheme. The re- 
sult follows naturally by combining the usual DFT based 
friction tensor j2 II] with the relationship between the fric- 
tion and fluctua ting forces in a non-adiabatic Newns- 
Anderson Model. 17] In principle, this scheme allows one 
to model non-adiabatic dynamics at metal surfaces by 
Langevin dynamics with ab-initio non-Markovian friction 
and fluctuating forces. 
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Appendix A: Path integral representation of the 
reduced density matrix 

In this appendix we give a path integral representation 
of the reduced density matrix Eq. in a coordinate 
basis. We will focus on the probabilistic interpretation 
of the path integral which leads to a Gaussian distributed 
classical Langevin equation. 

In a coordinate representation the reduced density ma- 
trix is 

Pred{x,y;t) = {x\Trei[p{t)]\y). (Al) 

and the diagonal elements give the probabilities of finding 
the adsorbate at a particular position regardless of the 
state of the electronic system. As shown in Ref. [13] 
the reduced density matrix of the Hamiltonian ([!]) can 
be represented as a double path integral: 



Predix,y;t) = J dxodyo{xo\po\yo) J P[a;(i')]2?[2/(i')]e^^="'"^''^'"^*'^'/', (A2) 
with the effective action given by 

Seff[x{t'),y{t')] = So[x{t')] - So[y{t')] - fdt' f dt"v{t'Ut' - t")u{t") (A3) 

Jo Jo 

+ \ I dt' [ dt"v{t')K{t' ^t")v{t"), 



2hJo 

where u{t) — x{t)/2 + y{t)/2, v{t) — x{t) — y{t) and -qlt) and K{t) are given by Eqs. (fTOj) and (fT3| respectively. 
With a quadratic potential the non-interacting action is given by So[x{t')] = M J^* dt'{x^ — uj'^x'^)/2. Changing to 
coordinates u and v and performing a partial integration on the kinetic term then gives for the diagonal part of the 
density matrix: 

Prea{u-t) = J duodpoV{uo,po) J P[u(^')]P[^.(^')]e-^^° '^'^^''^''^^'^-^^o (A4) 
where V{xo,Pq) is the Wigner distribution of pq, 

£,{t) = Mu{t) + Mulu^it) + [ dt'r]{t - t')u{t'), (A5) 

J to 

and u(t') have the additional constraint that uq = po/m. For a non-quadratic potential V{x), one is forced to make 
a semiclassical expansion of the potential to second order and the exponential in Eq. (jA4l) would contain additional 
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terms of order 0{v^V"'). Similarly, with a nonlinear interaction Hj in ([T]) one can perform a semiclassical expansion 
of the frictional terms which leads to a position dependence in Eq. ([T^. 

Without the quadratic term in v{t'), the density matrix ()A4|) would give a delta functional in £^{t) and the dynamics 
would be governed by a classical equation of motion with dynamical friction function r]{t). However the last term in 
the exponential of (|A4I) gives rise to a Gaussian broadening of the classical path. To see this explicitly we " complete" 
the square in the exponential and perform the path integral in v{t') which gives 



Prediu;t)oc / duodpoViuo,po) / V[u{t')]e-iIodt'dt"at')K-\t'-t")at")^ 



where K ^ solves 



dt"K-\t' - t")K{t" - t'") - 5{t' - t'"). 



(A6) 



(A7) 



The exponential in (|A6I) can be interpreted as the probability density of taking the path u{t') given the endpoints uq 
and u{t) and the initial velocity iio — po/m. It has a maximum at ^(t) = corresponding to the classical dynamics 
and the classical path is broadened by K~^. However, it will be more convenient to consider the probability density of 
^(t) which obviously has dimensions of a force. It is then necessary to change the path integral measure from T)[u{t')] 
to 23 [^(i')] and it can be shown that the Jacobian of this transformation is independent of The two-point 

correlation function of (^{t) can then be calculated by 



^^^^ ^^^^^ ^ J V[^{t')]e-i fo dt'dt"at')K-\t'-t")at")-f^ dt'j{t')i{t') 



.1=0 



!^dt'dt"i{f)K-^t'-t")£,{t") 



5J{ti)5J{t2) 



i dt'dt" J(t')K(t' -t")J{t") 



e2 Jo 



K{h-t2). 



(A8) 



.7=0 



This is the most compact way of specifying the statisti- 
cal properties of £^{t) and Eq. (jASP can be regarded as 
a classical equation motion with a stochastic Gaussian 
distributed force ^(t). 



Appendix B: Discretization of the correlation 
function 



mation can be obtained by a Cholesky decomposition of 
Cij such that = Y^j ^ijCj^ where J2j LijLkj = Cik- 

The stochastic force appearing in the Langevin equa- 
tion can be regarded as an infinite number of stochastic 
variables; one for each point in time from to t. Thus, 
to obtain an expression for the fiuctuation force in a time 
interval At, we need the statistical properties of the in- 
tegrals 



To sample a correlated "force path" £^{t) we need to 
discretize the correlation function and diagonalize the 
resulting correlation matrix. For a set of Gaussian dis- 
tributed stochastic variables {^i} with probability distri- 
bution 



(Bl) 



The correlation matrix Cij can be assumed symmetric 
without loss of generality. Hence, there exist a diagonal 
basis of uncorrelated variables {^^} which can be sampled 
from independent normalized Gaussians. The transfor- 



1 /•(i+l)At 



(B2) 



From the theory of multivariate Gaussian distributions it 
is readily shown that the set of these integrals are Gaus- 
sian distributed with the correlation matrix: 



-, /.(i+l)At /-(j + ljAt 

C,, = -r72 dt' dt"K{t' -t"), (B3) 

JiAt JiAt 



and this is the expression used when calculating molecu- 
lar trajectories using non-Markovian Langevin dynamics. 
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